function brt2rad_f(np, wavenumbers, brt) result(rad)
    use iasing_osse
    implicit none
    ! Arguments
    integer(JPI4), intent(in)   :: np
    real(JPR4), dimension(np), intent(in)   :: wavenumbers, brt
    real(JPR4), dimension(np)   :: rad

    ! Auxiliar variables
    real(JPR8), dimension(np)   :: rad8, brt8
    real(JPR8), dimension(np)   :: wnu8_scaled
    real(JPR8), parameter       :: scale = 100

    wnu8_scaled = real(wavenumbers, 8)*scale
    brt8        = real(brt, 8)
    
    rad8 =  (scale * CST_2HC2 * wnu8_scaled**3) / &
            (exp(CST_HCK*wnu8_scaled/brt8)-1)

    rad  = real(rad8, 4)

end function brt2rad_f
